<?php

// 定义 sigmoid 函数
function sigmoid($z)
{
    return 1 / (1 + exp(-1 * $z));
}

// 以 v 填充二维数组 r x c
function fill2d($r, $c, $v)
{
    $arr = [];
    for ($i = 0; $i < $r; $i++) {
        $arr[] = array_pad([], $c, $v);
    }
    return $arr;
}

// 矩阵按行转向量
function mat2vec(array $mat)
{
    $m = count($mat);
    $n = count($mat[0]);
    $vec = [];
    for ($i = 0; $i < $m; $i++) {
        for ($j = 0; $j < $n; $j++) {
            $vec[] = $mat[$i][$j];
        }
    }
    return $vec;
}

// 向量按行转矩阵 r x c
function vec2mat(array $vec, int $r, int $c)
{
    if (count($vec) != $r * $c) {
        throw new Exception("Number of vector cell don't equal r x c");
    }
    $mat = [];
    for ($i = 0; $i < $r; $i++) {
        $row = [];
        for ($j = 0; $j < $c; $j++) {
            $row[] = $vec[$i * $c + $j];
        }
        $mat[] = $row;
    }
    return $mat;
}

// ##### 前向传播 ################
function fp(array $a1, array $w)
{
    $m = count($a1);
    $n = count($a1[0]);

    $ws = array_chunk($w, 6);
    $w1 = vec2mat($ws[0], 2, 3);
    $w2 = vec2mat($ws[1], 2, 3);

    $A1 = fill2d($m, 3, 1);
    $a2 = fill2d($m, 2, 0);
    $A2 = fill2d($m, 3, 1);
    $a3 = fill2d($m, 2, 0);

    // 输入加偏执单元
    for ($i = 0; $i < $m; $i++) {
        $A1[$i][0] = 1;
        $A1[$i][1] = $a1[$i][0];
        $A1[$i][2] = $a1[$i][1];
    }

    // compute A2, a3
    for ($i = 0; $i < $m; $i++) {
        // 1 -> 2 层
        // 加权和
        $z11 = $A1[$i][0] * $w1[0][0] + $A1[$i][1] * $w1[0][1] + $A1[$i][2] * $w1[0][2];
        $z12 = $A1[$i][0] * $w1[1][0] + $A1[$i][1] * $w1[1][1] + $A1[$i][2] * $w1[1][2];

        $a2[$i][0] = sigmoid($z11);
        $a2[$i][1] = sigmoid($z12);

        // A2 加入偏置单元
        $A2[$i][0] = 1;
        $A2[$i][1] = $a2[$i][0];
        $A2[$i][2] = $a2[$i][1];

        // 2 -> 3 层
        $z21 = $A2[$i][0] * $w2[0][0] + $A2[$i][1] * $w2[0][1] + $A2[$i][2] * $w2[0][2];
        $z22 = $A2[$i][0] * $w2[1][0] + $A2[$i][1] * $w2[1][1] + $A2[$i][2] * $w2[1][2];
        $a3[$i][0] = sigmoid($z21);
        $a3[$i][1] = sigmoid($z22);
    }

    return $a3;
}

/**
 * 这不是一个通用函数，它只能处理:
 *     输入特征数为 2
 *     一个隐藏层，特征数为 2
 *     输出层特征数为 2
 * 的数据
 */
function compute(array $a1, array $y, array $w, $lambda = 0)
{
    $m = count($a1);
    $n = count($a1[0]);

    // 参数向量转矩阵
    $ws = array_chunk($w, 6);
    $w1 = vec2mat($ws[0], 2, 3);
    $w2 = vec2mat($ws[1], 2, 3);

    $A1 = fill2d($m, 3, 1);
    $a2 = fill2d($m, 2, 0);
    $A2 = fill2d($m, 3, 1);
    $a3 = fill2d($m, 2, 0);

    // 输入加偏执单元
    for ($i = 0; $i < $m; $i++) {
        $A1[$i][0] = 1;
        $A1[$i][1] = $a1[$i][0];
        $A1[$i][2] = $a1[$i][1];
    }

    // compute A2, a3
    for ($i = 0; $i < $m; $i++) {
        // 1 -> 2 层
        // 加权和
        $z11 = $A1[$i][0] * $w1[0][0] + $A1[$i][1] * $w1[0][1] + $A1[$i][2] * $w1[0][2];
        $z12 = $A1[$i][0] * $w1[1][0] + $A1[$i][1] * $w1[1][1] + $A1[$i][2] * $w1[1][2];

        $a2[$i][0] = sigmoid($z11);
        $a2[$i][1] = sigmoid($z12);

        // A2 加入偏置单元
        $A2[$i][0] = 1;
        $A2[$i][1] = $a2[$i][0];
        $A2[$i][2] = $a2[$i][1];

        // 2 -> 3 层
        $z21 = $A2[$i][0] * $w2[0][0] + $A2[$i][1] * $w2[0][1] + $A2[$i][2] * $w2[0][2];
        $z22 = $A2[$i][0] * $w2[1][0] + $A2[$i][1] * $w2[1][1] + $A2[$i][2] * $w2[1][2];
        $a3[$i][0] = sigmoid($z21);
        $a3[$i][1] = sigmoid($z22);
    }

    // 计算损失
    $J = 0;
    for ($i = 0; $i < $m; $i++) {
        $J += pow($a3[$i][0] - $y[$i][0], 2) + pow($a3[$i][1] - $y[$i][1], 2);
    }
    $J = $J / (2 * $m);

    // 加入正则化项
    $r = 0;
    for ($i = 0; $i < 2; $i++) {
        // $j 从1开始也可以
        for ($j = 0; $j < 3; $j++) {
            $r += pow($w1[$i][$j], 2);
            $r += pow($w2[$i][$j], 2);
        }
    }
    $r = $r * $lambda / (2 * $m);

    // 结构风险最小化损失
    $Jsrm = $J + $r;

    // ###### 后向传播 ################
    // 初始化 J 对 w2 的偏导
    $p2 = fill2d(2, 3, 0);
    // 初始化 J 对 w1 的偏导
    $p1 = fill2d(2, 3, 0);

    for ($i = 0; $i < $m; $i++) {
        $delta1 = 1 / $m * ($a3[$i][0] - $y[$i][0]) * $a3[$i][0] * (1 - $a3[$i][0]);
        $delta2 = 1 / $m * ($a3[$i][1] - $y[$i][1]) * $a3[$i][1] * (1 - $a3[$i][1]);
        $p2[0][0] += $delta1 * $A2[$i][0];
        $p2[0][1] += $delta1 * $A2[$i][1];
        $p2[0][2] += $delta1 * $A2[$i][2];
        $p2[1][0] += $delta2 * $A2[$i][0];
        $p2[1][1] += $delta2 * $A2[$i][1];
        $p2[1][2] += $delta2 * $A2[$i][2];
    }
    // 加入正则化项的偏导
    $p2[0][0] += $lambda / $m * $p2[0][0]; // 这一项不加也可以
    $p2[0][1] += $lambda / $m * $p2[0][1];
    $p2[0][2] += $lambda / $m * $p2[0][2];
    $p2[1][0] += $lambda / $m * $p2[1][0]; // 这一项不加也可以
    $p2[1][1] += $lambda / $m * $p2[1][1];
    $p2[1][2] += $lambda / $m * $p2[1][2];

    for ($i = 0; $i < $m; $i++) {
        $delta1 = 1 / $m * (
            ($a3[$i][0] - $y[$i][0]) * $a3[$i][0] * (1 - $a3[$i][0]) * $w2[0][1] +
            ($a3[$i][1] - $y[$i][1]) * $a3[$i][1] * (1 - $a3[$i][1]) * $w2[1][1]
        ) * $A2[$i][1] * (1 - $A2[$i][1]);

        $delta2 = 1 / $m * (
            ($a3[$i][0] - $y[$i][0]) * $a3[$i][0] * (1 - $a3[$i][0]) * $w2[0][2] +
            ($a3[$i][1] - $y[$i][1]) * $a3[$i][1] * (1 - $a3[$i][1]) * $w2[1][2]
        ) * $A2[$i][2] * (1 - $A2[$i][2]);

        $p1[0][0] += $delta1 * $A1[$i][0];
        $p1[0][1] += $delta1 * $A1[$i][1];
        $p1[0][2] += $delta1 * $A1[$i][2];
        $p1[1][0] += $delta2 * $A1[$i][0];
        $p1[1][1] += $delta2 * $A1[$i][1];
        $p1[1][2] += $delta2 * $A1[$i][2];
    }

    // 加入正则化项的偏导
    $p1[0][0] += $lambda / $m * $p1[0][0]; // 这一项不加也可以
    $p1[0][1] += $lambda / $m * $p1[0][1];
    $p1[0][2] += $lambda / $m * $p1[0][2];
    $p1[1][0] += $lambda / $m * $p1[1][0]; // 这一项不加也可以
    $p1[1][1] += $lambda / $m * $p1[1][1];
    $p1[1][2] += $lambda / $m * $p1[1][2];

    // 偏导矩阵转向量
    $p = array_merge(mat2vec($p1), mat2vec($p2));

    return ["J" => $J, "P" => $p, "HatY" => $a3];
}

// ############# 梯度下降求权重 ##############
    // w <- w - alpha * (partial J) / (partial w)
function gradDesc(array $w, float $alpha, array $partial)
{
    $n = count($w);
    for ($i = 0; $i < $n; $i++) {
            $w[$i] = $w[$i] - $alpha * $partial[$i];
    }
    return $w;
}
